NOTE: from level 4/5 onwards I will execute most steps locally, as the datasets are not as large anymore.

1 Introduction

Here, we will include the annotation of the Cytotoxic cells (level 4).

1.1 Load packages

library(Seurat)
library(SeuratWrappers)
library(harmony)
library(clusterProfiler)
library(org.Hs.eg.db)
library(tidyverse)

1.2 Parameters

# Paths
path_to_obj <- "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/results/R_objects/level_3/Cytotoxic/Cytotoxic_clustered_level_3.rds"
path_to_naive_cd8 <- "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/results/R_objects/level_4/CD4_T/naive_CD8_T_subsetted_level_4.rds"
path_to_save <- "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/results/R_objects/level_4/Cytotoxic/Cytotoxic_annotated_level_4.rds"


# Colors
color_palette <-  c("#1CFFCE", "#90AD1C", "#C075A6", "#85660D", "#5A5156", "#AA0DFE",   
                    "#F8A19F", "#F7E1A0", "#1C8356", "#FEAF16", "#822E1C", "#C4451C",   
                    "#1CBE4F", "#325A9B", "#F6222E", "#FE00FA", "#FBE426", "#16FF32", 
                    "black",   "#3283FE", "#B00068", "#DEA0FD", "#B10DA1", "#E4E1E3",   
                    "#90AD1C", "#FE00FA", "#85660D", "#3B00FB", "#822E1C", "coral2", 
                    "#1CFFCE", "#1CBE4F", "#3283FE", "#FBE426", "#F7E1A0", "#325A9B",   
                    "#2ED9FF", "#B5EFB5", "#5A5156", "#DEA0FD", "#FEAF16", "#683B79",   
                    "#B10DA1", "#1C7F93", "#F8A19F", "dark orange", "#FEAF16", "#FBE426",  
                    "Brown")

# Functions
source("~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/bin/utils.R")

1.3 Load data

seurat <- readRDS(path_to_obj)
seurat

## An object of class Seurat 
## 37378 features across 10998 samples within 1 assay 
## Active assay: RNA (37378 features, 0 variable features)
##  3 dimensional reductions calculated: pca, umap, harmony

DimPlot(seurat, cols = color_palette)

naive_cd8 <- readRDS(path_to_naive_cd8)

2 Important literature

3 Key steps

Isolate CD8 T cells, Including the naive CD8 T cells that clustered with CD4 T. This will help us differentiate between LEF1+SELL+IL7R+ from LEF1+SELL+IL7R-. Subcluster 13 to find CD8 T effector cytotoxic and CD16+ NK. Watch out with doublets with CD4 T.
Subcluster TIM3 DN into two clusters.
Subset NK and ILC and annotate.
Differential expression analysis (6 vs 11, 2 vs 6, 2 vs 11).
Subcluster cluster 7 (exhausted) into 2.
Subcluster cluster 5 (MAIT) into 2.

4 Preliminar annotation

https://www.nature.com/articles/s41467-018-06318-7 s * 0: * 1: CD8 T effector memory * 2: CD8 T naive-like (Li H) OR CD8 T memory (Sade-Feldman) OR CD8 T stem cell-like (Gattinoni) * 3: CD8 T stem cell-like? (PTPRC) * 4: TIM3 DN * 5: MAIT cells https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3192229/. CD69+ –> tissue-resident MAIT cells? * 6: gamma-delta T cells effector? * 7: CD8 T exhausted * 8: ILC * 9: CD16- NK
* 10: gamma-delta T cells * 11: TSHZ2+CD130+_doublets * 12: CD226 effector? * 13: CD8 T cytotoxic effector + CD16+ NK * 14: interferon * 15: Doublets * 16: poor-quality

Resident cells? Central memory?

5 Subcluster

5.1 Cluster 4 (TIM3 DN)

Idents(seurat) <- "seurat_clusters"
seurat <- FindSubCluster(seurat, "4", "RNA_snn", "TIM3", resolution = 0.2)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 958
## Number of edges: 39455
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8138
## Number of communities: 2
## Elapsed time: 0 seconds

Idents(seurat) <- "TIM3"
DimPlot(seurat, group.by = "TIM3", cols = color_palette)

markers_tim3 <- find_markers_subclusters(seurat, "TIM3", "^4_")
DT::datatable(markers_tim3$`4_0`)

DT::datatable(markers_tim3$`4_1`)

seurat$annotation_level_3 <- seurat$TIM3
seurat$annotation_level_3[seurat$annotation_level_3 == "4_0"] <- "4_CD200-TIM3+ DN"
seurat$annotation_level_3[seurat$annotation_level_3 == "4_1"] <- "4_CD200+TIM3+ DN"

5.2 Cluster 5 (MAIT cells)

Idents(seurat) <- "annotation_level_3"
seurat <- FindSubCluster(seurat, "5", "RNA_snn", "MAIT", resolution = 0.2)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 841
## Number of edges: 24088
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8280
## Number of communities: 3
## Elapsed time: 0 seconds

Idents(seurat) <- "MAIT"
DimPlot(seurat, group.by = "MAIT", cols = color_palette)

markers_mait <- find_markers_subclusters(seurat, "MAIT", "^5_")
DT::datatable(markers_mait$`5_0`)

DT::datatable(markers_mait$`5_1`)

DT::datatable(markers_mait$`5_2`)

seurat$annotation_level_3 <- seurat$MAIT
seurat$annotation_level_3[seurat$annotation_level_3 == "5_0"] <- "5_MAIT1_NK_KLRF1"
seurat$annotation_level_3[seurat$annotation_level_3 == "5_1"] <- "5_MAIT2_NK_KLRF1"
seurat$annotation_level_3[seurat$annotation_level_3 == "5_2"] <- "5_MAIT3"
DimPlot(seurat, group.by = "annotation_level_3", cols = color_palette)

5.3 Cluster 6 (effector gamma-delta T cells)

Idents(seurat) <- "annotation_level_3"
seurat <- FindSubCluster(seurat, "6", "RNA_snn", "gamma_delta", resolution = 0.2)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 639
## Number of edges: 25223
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8107
## Number of communities: 3
## Elapsed time: 0 seconds

Idents(seurat) <- "gamma_delta"
DimPlot(seurat, group.by = "gamma_delta", cols = color_palette)

markers_gamma_delta <- find_markers_subclusters(seurat, "gamma_delta", "^6_")
DT::datatable(markers_gamma_delta$`6_0`)

DT::datatable(markers_gamma_delta$`6_1`)

DT::datatable(markers_gamma_delta$`6_2`)

seurat$annotation_level_3 <- seurat$gamma_delta
seurat$annotation_level_3[seurat$annotation_level_3 == "6_0"] <- "6_effector_gamma_delta/NKT 1"
seurat$annotation_level_3[seurat$annotation_level_3 == "6_1"] <- "6_effector_gamma_delta/NKT 2"
seurat$annotation_level_3[seurat$annotation_level_3 == "6_2"] <- "6_effector_gamma_delta/NKT 3"
DimPlot(seurat, group.by = "annotation_level_3", cols = color_palette)

5.4 Cluster 7 (exhausted + doublets)

Idents(seurat) <- "annotation_level_3"
seurat <- FindSubCluster(seurat, "7", "RNA_snn", "exhausted_and_doublets", resolution = 0.1)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 574
## Number of edges: 19639
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9220
## Number of communities: 2
## Elapsed time: 0 seconds

Idents(seurat) <- "exhausted_and_doublets"
DimPlot(seurat, group.by = "exhausted_and_doublets", cols = color_palette)

markers_exh_doubl <- find_markers_subclusters(seurat, "exhausted_and_doublets", "^7_")
DT::datatable(markers_exh_doubl$`7_0`)

DT::datatable(markers_exh_doubl$`7_1`)

seurat$annotation_level_3 <- seurat$exhausted_and_doublets
seurat$annotation_level_3[seurat$annotation_level_3 == "7_0"] <- "7_CD8 T exhausted"
seurat$annotation_level_3[seurat$annotation_level_3 == "7_1"] <- "7_doublets"
DimPlot(seurat, group.by = "annotation_level_3", cols = color_palette)

5.5 Cluster 8 (ILC)

Idents(seurat) <- "annotation_level_3"
seurat <- FindSubCluster(seurat, "8", "RNA_snn", "ilc", resolution = 0.6)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 514
## Number of edges: 23484
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6772
## Number of communities: 3
## Elapsed time: 0 seconds

Idents(seurat) <- "ilc"
DimPlot(seurat, group.by = "ilc", cols = color_palette)

markers_ilc <- find_markers_subclusters(seurat, "ilc", "^8_")
DT::datatable(markers_ilc$`8_0`)

DT::datatable(markers_ilc$`8_1`)

DT::datatable(markers_ilc$`8_2`)

seurat$annotation_level_3 <- seurat$ilc
seurat$annotation_level_3[seurat$annotation_level_3 == "8_0"] <- "8_ILC2_RORA"
seurat$annotation_level_3[seurat$annotation_level_3 == "8_1"] <- "8_ILC1_CD200R1"
seurat$annotation_level_3[seurat$annotation_level_3 == "8_2"] <- "8_ILC3_RORC"
DimPlot(seurat, group.by = "annotation_level_3", cols = color_palette)

5.6 Cluster 9 (NK)

Idents(seurat) <- "annotation_level_3"
seurat <- FindSubCluster(seurat, "9", "RNA_snn", "nk", resolution = 0.6)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 411
## Number of edges: 19675
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6208
## Number of communities: 2
## Elapsed time: 0 seconds

Idents(seurat) <- "nk"
DimPlot(seurat, group.by = "nk", cols = color_palette)

markers_nk <- find_markers_subclusters(seurat, "nk", "^9_")
DT::datatable(markers_nk$`9_0`)

DT::datatable(markers_nk$`9_1`)

seurat$annotation_level_3 <- seurat$nk
seurat$annotation_level_3[seurat$annotation_level_3 == "9_0"] <- "9_CD16+CD56- NK"
seurat$annotation_level_3[seurat$annotation_level_3 == "9_1"] <- "9_CD16+CD56+ NK"
DimPlot(seurat, group.by = "annotation_level_3", cols = color_palette)

5.7 Cluster 13 (CTL + CD16+ NK)

Idents(seurat) <- "annotation_level_3"
seurat <- FindSubCluster(seurat, "13", "RNA_snn", "ctl_cd16_nk", resolution = 0.3)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 142
## Number of edges: 3221
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7630
## Number of communities: 2
## Elapsed time: 0 seconds

Idents(seurat) <- "ctl_cd16_nk"
DimPlot(seurat, group.by = "ctl_cd16_nk", cols = color_palette)

markers_ctl_cd16_nk <- find_markers_subclusters(seurat, "ctl_cd16_nk", "^13_")
DT::datatable(markers_ctl_cd16_nk$`13_0`)

DT::datatable(markers_ctl_cd16_nk$`13_1`)

seurat$annotation_level_3 <- seurat$ctl_cd16_nk
seurat$annotation_level_3[seurat$annotation_level_3 == "13_0"] <- "13_CTL"
seurat$annotation_level_3[seurat$annotation_level_3 == "13_1"] <- "13_CD16_NK"
DimPlot(seurat, group.by = "annotation_level_3", cols = color_palette)

5.8 Finish annotation at this level

seurat$annotation_level_3[seurat$annotation_level_3 == "0"] <- "0_TBD"
seurat$annotation_level_3[seurat$annotation_level_3 == "1"] <- "1_CD8_T_effector_memory"
seurat$annotation_level_3[seurat$annotation_level_3 == "2"] <- "2_naive_like/memory/stem-like"
seurat$annotation_level_3[seurat$annotation_level_3 == "3"] <- "3_naive_like/memory/stem-like"
seurat$annotation_level_3[seurat$annotation_level_3 == "10"] <- "10_gamma_delta_T"
seurat$annotation_level_3[seurat$annotation_level_3 == "11"] <- "11_TSHZ2_CD130_doublets"
seurat$annotation_level_3[seurat$annotation_level_3 == "12"] <- "12_CD226_CD55"
seurat$annotation_level_3[seurat$annotation_level_3 == "14"] <- "14_IFNG_CD8_T"
seurat$annotation_level_3[seurat$annotation_level_3 == "15"] <- "15_doublets"
seurat$annotation_level_3[seurat$annotation_level_3 == "16"] <- "16_poor_quality"
Idents(seurat) <- "annotation_level_3"
DimPlot(seurat, group.by = "annotation_level_3", pt.size = 0.75, cols = color_palette)

6 Exclude cells

We will exclude:

Cluster 15: doublets
Cluster 16: poor-quality
A subset of cluster 7 which are doublets with CD4
Cells from King et al. or multiome+old adult (decided jointly with other omics).

donors_to_exclude <- str_subset(unique(seurat$donor_id) , "BCP")
to_exclude <- 
  (seurat$annotation_level_3 %in% c("15_doublets", "16_poor_quality", "7_doublets")) |
  (seurat$assay == "5P") |
  (seurat$assay == "multiome" & seurat$donor_id == "BCLL-2-T") |
  (seurat$donor_id %in% donors_to_exclude)
seurat <- subset(seurat, cells = colnames(seurat)[!to_exclude])
seurat

## An object of class Seurat 
## 37378 features across 8398 samples within 1 assay 
## Active assay: RNA (37378 features, 0 variable features)
##  3 dimensional reductions calculated: pca, umap, harmony

DimPlot(seurat, cols = color_palette, group.by = "annotation_level_3")

7 Create datasets level 4

7.1 CD8 T cells

seurat <- subset(seurat, annotation_level_3 != "0_TBD")
cd8_levels <- c(
  "1_CD8_T_effector_memory",
  "10_gamma_delta_T",
  "11_TSHZ2_CD130_doublets",
  "12_CD226_CD55",
  "13_CTL",
  "14_IFNG_CD8_T",
  "2_naive_like/memory/stem-like",
  "3_naive_like/memory/stem-like",
  "5_MAIT1_NK_KLRF1",
  "5_MAIT2_NK_KLRF1",
  "5_MAIT3",
  "6_effector_gamma_delta/NKT 1",
  "6_effector_gamma_delta/NKT 2",
  "6_effector_gamma_delta/NKT 3",
  "7_CD8 T exhausted"
)
Idents(seurat) <- "annotation_level_3"
cd8 <- subset(seurat, idents = cd8_levels)


# Reprocess
cd8_list <- SplitObject(cd8, split.by = "assay")
cd8_list <- cd8_list[c("3P", "multiome")]
cd8_list <- purrr::map(
  cd8_list,
  FindVariableFeatures,
  nfeatures = 5000
)
hvg <- purrr::map(cd8_list, VariableFeatures)
shared_hvg <- intersect(hvg$`3P`, hvg$multiome)
cd8 <- cd8 %>%
  ScaleData(features = shared_hvg) %>%
  RunPCA(features = shared_hvg) %>%
  RunHarmony(group.by.vars = "assay", reduction ="pca", dims = 1:25) %>%
  RunUMAP(reduction = "harmony", dims = 1:25)
DimPlot(cd8, group.by = "annotation_level_3", cols = color_palette)

# Recluster
cd8 <- FindNeighbors(cd8, reduction = "harmony", dims = 1:25)
cd8 <- FindClusters(cd8, resolution = 1.1)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5272
## Number of edges: 222040
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8086
## Number of communities: 15
## Elapsed time: 0 seconds

DimPlot(cd8, cols = color_palette, pt.size = 0.75)

# Visualize proportions
proportions_cd8_df <- find_proportions_df(
  cd8,
  x = "seurat_clusters",
  fill = "annotation_level_3"
)
proportions_cd8_gg <- plot_stacked_barplot(
  proportions_cd8_df,
  x = "seurat_clusters",
  fill = "annotation_level_3",
  colors = color_palette
)
proportions_cd8_gg

# Markers
markers_cd8 <- FindAllMarkers(
  cd8,
  only.pos = TRUE,
  logfc.threshold = 0.6
)
markers_cd8 <- markers_cd8 %>%
  group_by(cluster) %>%
  arrange(desc(avg_log2FC), .by_group = TRUE)
DT::datatable(markers_cd8)

Clusters 9, 10 and 11 are doublets and/or poor-quality cells:

selected_cells <- colnames(cd8)[!(cd8$seurat_clusters %in% c("9", "10", "11"))]
cd8 <- subset(cd8, cells = selected_cells)
DimPlot(cd8, cols = color_palette)

# Reprocess
cd8_list <- SplitObject(cd8, split.by = "assay")
cd8_list <- cd8_list[c("3P", "multiome")]
cd8_list <- purrr::map(
  cd8_list,
  FindVariableFeatures,
  nfeatures = 5000
)
hvg <- purrr::map(cd8_list, VariableFeatures)
shared_hvg <- intersect(hvg$`3P`, hvg$multiome)
cd8 <- cd8 %>%
  ScaleData(features = shared_hvg) %>%
  RunPCA(features = shared_hvg) %>%
  RunHarmony(group.by.vars = "assay", reduction ="pca", dims = 1:25) %>%
  RunUMAP(reduction = "harmony", dims = 1:25)
DimPlot(cd8, group.by = "annotation_level_3", cols = color_palette)

# Recluster
cd8 <- FindNeighbors(cd8, reduction = "harmony", dims = 1:25)
cd8 <- FindClusters(cd8, resolution = 1)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 4626
## Number of edges: 198873
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8230
## Number of communities: 15
## Elapsed time: 0 seconds

DimPlot(cd8, cols = color_palette, pt.size = 0.75)

# Visualize proportions
proportions_cd8_df <- find_proportions_df(
  cd8,
  x = "seurat_clusters",
  fill = "annotation_level_3"
)
proportions_cd8_gg <- plot_stacked_barplot(
  proportions_cd8_df,
  x = "seurat_clusters",
  fill = "annotation_level_3",
  colors = color_palette
)+ ggtitle ("")
proportions_cd8_gg

# Markers
markers_cd8 <- FindAllMarkers(
  cd8,
  only.pos = TRUE,
  logfc.threshold = 0.6
)
markers_cd8 <- markers_cd8 %>%
  group_by(cluster) %>%
  arrange(desc(avg_log2FC), .by_group = TRUE)
DT::datatable(markers_cd8)

markers_cd8_list <- purrr::map(unique(markers_cd8$cluster), function(x) {
  df <- markers_cd8[markers_cd8$cluster == x, ]
  df
})
names(markers_cd8_list) <- unique(markers_cd8$cluster)

7.2 ILC/NK

ilc_nk_levels <- c(
  "8_ILC1_CD200R1",
  "8_ILC2_RORA",
  "8_ILC3_RORC",
  "9_CD16+CD56+ NK",
  "9_CD16+CD56- NK",
  "13_CD16_NK"
)
Idents(seurat) <- "annotation_level_3"
ilc_nk <- subset(seurat, idents = ilc_nk_levels)
DimPlot(ilc_nk)

ilc_nk_list <- SplitObject(ilc_nk, split.by = "assay")
ilc_nk_list <- ilc_nk_list[c("3P", "multiome")]
ilc_nk_list <- purrr::map(
  ilc_nk_list,
  FindVariableFeatures,
  nfeatures = 5000
)
hvg <- purrr::map(ilc_nk_list, VariableFeatures)
shared_hvg <- intersect(hvg$`3P`, hvg$multiome)
ilc_nk <- ilc_nk %>%
  ScaleData(features = shared_hvg) %>%
  RunPCA(features = shared_hvg) %>%
  RunHarmony(group.by.vars = "assay", reduction ="pca", dims = 1:20) %>%
  RunUMAP(reduction = "harmony", dims = 1:20)
DimPlot(ilc_nk, group.by = "annotation_level_3", cols = color_palette)

ilc_nk <- FindNeighbors(ilc_nk, dims = 1:20, reduction = "harmony")
ilc_nk <- FindClusters(ilc_nk, resolution = 0.6)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 719
## Number of edges: 27685
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8145
## Number of communities: 7
## Elapsed time: 0 seconds

DimPlot(ilc_nk)

# Visualize proportions
proportions_ilc_nk_df <- find_proportions_df(
  ilc_nk,
  x = "seurat_clusters",
  fill = "annotation_level_3"
)
proportions_ilc_nk_gg <- plot_stacked_barplot(
  proportions_ilc_nk_df,
  x = "seurat_clusters",
  fill = "annotation_level_3",
  colors = color_palette
)
proportions_ilc_nk_gg

FeaturePlot(ilc_nk, "CD3D")

We can subset cluster 5 (doublets with T cells):

ilc_nk <- subset(ilc_nk, seurat_clusters != "5")
DimPlot(ilc_nk, cols = color_palette, group.by = "seurat_clusters")

ilc_nk <- FindNeighbors(ilc_nk, reduction = "harmony", dims = 1:20)
Idents(ilc_nk) <- "seurat_clusters"
ilc_nk <- FindSubCluster(ilc_nk, "1", "RNA_snn", "ilc_nk_1", resolution = 0.2)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 149
## Number of edges: 5024
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8096
## Number of communities: 2
## Elapsed time: 0 seconds

Idents(ilc_nk) <- "ilc_nk_1"
DimPlot(ilc_nk, group.by = "ilc_nk_1", cols = color_palette)

ilc_nk$seurat_clusters <- factor(
  ilc_nk$ilc_nk_1,
  levels = c("0", "1_0", "1_1", "2", "3", "4", "6")
)
Idents(ilc_nk) <- "seurat_clusters"


# Visualize proportions
proportions_ilc_nk_df <- find_proportions_df(
  ilc_nk,
  x = "seurat_clusters",
  fill = "annotation_level_3"
)
proportions_ilc_nk_gg <- plot_stacked_barplot(
  proportions_ilc_nk_df,
  x = "seurat_clusters",
  fill = "annotation_level_3",
  colors = color_palette
) + ggtitle ("")
proportions_ilc_nk_gg

# Markers
markers_ilc_nk <- FindAllMarkers(
  ilc_nk,
  only.pos = TRUE,
  logfc.threshold = 0.6
)
markers_ilc_nk <- markers_ilc_nk %>%
  group_by(cluster) %>%
  arrange(desc(avg_log2FC), .by_group = TRUE)
DT::datatable(markers_ilc_nk)

markers_ilc_nk_list <- purrr::map(unique(markers_ilc_nk$cluster), function(x) {
  df <- markers_ilc_nk[markers_ilc_nk$cluster == x, ]
  df
})
names(markers_ilc_nk_list) <- unique(markers_ilc_nk$cluster)

8 Save

# CD8 T
input_shiny_cd8 <- seurat2shiny(cd8, slot = "data", reduction = "umap")
saveRDS(
  input_shiny_cd8$expression,
  "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/results/R_objects/level_4/Cytotoxic/CD8_T/CD8_T_expression_to_shiny_app_level_4.rds"
)
saveRDS(
  input_shiny_cd8$metadata,
  "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/results/R_objects/level_4/Cytotoxic/CD8_T/CD8_T_metadata_to_shiny_app_level_4.rds"
)
saveRDS(cd8, "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/results/R_objects/level_4/Cytotoxic/CD8_T/CD8_T_clustered_level_4.rds")
openxlsx::write.xlsx(
  markers_cd8_list,
  "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/3-clustering/4-level_4/Cytotoxic/tmp/CD8_T_markers.xlsx"
)
umap_level_4_cd8 <- DimPlot(cd8, cols = color_palette, pt.size = 1)
ggsave(
  filename = "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/3-clustering/4-level_4/Cytotoxic/tmp/CD8_T_umap_level_4.png",
  plot = umap_level_4_cd8,
  width = 14,
  height = 12,
  units = "cm"
)
ggsave(
  filename = "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/3-clustering/4-level_4/Cytotoxic/tmp/CD8_T_stacked_barplot_preliminar_annotation_level_4.png",
  plot = proportions_cd8_gg,
  width = 16,
  height = 12,
  units = "cm"
)

# ILC/NK
input_shiny_ilc_nk <- seurat2shiny(ilc_nk, slot = "data", reduction = "umap")
saveRDS(
  input_shiny_ilc_nk$expression,
  "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/results/R_objects/level_4/Cytotoxic/ILC_NK/ILC_NK_expression_to_shiny_app_level_4.rds"
)
saveRDS(
  input_shiny_ilc_nk$metadata,
  "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/results/R_objects/level_4/Cytotoxic/ILC_NK/ILC_NK_metadata_to_shiny_app_level_4.rds"
)
saveRDS(ilc_nk, "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/results/R_objects/level_4/Cytotoxic/ILC_NK/ILC_NK_clustered_level_4.rds")
openxlsx::write.xlsx(
  markers_ilc_nk_list,
  "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/3-clustering/4-level_4/Cytotoxic/tmp/ILC_NK_markers.xlsx"
)
umap_level_4_ilc_nk <- DimPlot(ilc_nk, cols = color_palette, pt.size = 1)
ggsave(
  filename = "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/3-clustering/4-level_4/Cytotoxic/tmp/ILC_NK_umap_level_4.png",
  plot = umap_level_4_ilc_nk,
  width = 14,
  height = 12,
  units = "cm"
)
ggsave(
  filename = "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/3-clustering/4-level_4/Cytotoxic/tmp/ILC_NK_stacked_barplot_preliminar_annotation_level_4.png",
  plot = proportions_ilc_nk_gg,
  width = 16,
  height = 12,
  units = "cm"
)

# TIM3 DN
tim3_dn <- subset(seurat, idents = c("4_CD200-TIM3+ DN", "4_CD200+TIM3+ DN"))
saveRDS(tim3_dn, "~/Desktop/CNAG/mnt_clust/tonsil_atlas/current/scRNA-seq/results/R_objects/level_4/Cytotoxic/TIM3_DN/TIM3_DN_clustered_level_4.rds")

9 Session Information

sessionInfo()

## R version 4.0.4 (2021-02-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=es_ES.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=es_ES.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=es_ES.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] forcats_0.5.1          stringr_1.4.0          dplyr_1.0.6            purrr_0.3.4            readr_1.4.0            tidyr_1.1.3            tibble_3.1.2           ggplot2_3.3.3          tidyverse_1.3.1        org.Hs.eg.db_3.12.0    AnnotationDbi_1.52.0   IRanges_2.24.1         S4Vectors_0.28.1       Biobase_2.50.0         BiocGenerics_0.36.1    clusterProfiler_3.18.1 harmony_1.0            Rcpp_1.0.6             SeuratWrappers_0.3.0   SeuratObject_4.0.2     Seurat_4.0.3           BiocStyle_2.18.1      
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.1            reticulate_1.20       tidyselect_1.1.1      RSQLite_2.2.7         htmlwidgets_1.5.3     grid_4.0.4            BiocParallel_1.24.1   Rtsne_0.15            scatterpie_0.1.6      munsell_0.5.0         codetools_0.2-18      ica_1.0-2             DT_0.18               future_1.21.0         miniUI_0.1.1.1        withr_2.4.2           colorspace_2.0-1      GOSemSim_2.16.1       highr_0.9             knitr_1.33            rstudioapi_0.13       ROCR_1.0-11           tensor_1.5            DOSE_3.16.0           listenv_0.8.0         labeling_0.4.2        polyclip_1.10-0       bit64_4.0.5           farver_2.1.0          downloader_0.4        parallelly_1.26.0     vctrs_0.3.8           generics_0.1.0        xfun_0.23             R6_2.5.0              graphlayouts_0.7.1    rsvd_1.0.5            spatstat.utils_2.2-0  cachem_1.0.5          fgsea_1.16.0          assertthat_0.2.1      promises_1.2.0.1      scales_1.1.1          ggraph_2.0.5          enrichplot_1.10.2     gtable_0.3.0          globals_0.14.0        goftest_1.2-2         tidygraph_1.2.0       rlang_0.4.11          splines_4.0.4         lazyeval_0.2.2        spatstat.geom_2.1-0   broom_0.7.7          
##  [55] BiocManager_1.30.15   yaml_2.2.1            reshape2_1.4.4        abind_1.4-5           modelr_0.1.8          crosstalk_1.1.1       backports_1.2.1       httpuv_1.6.1          qvalue_2.22.0         tools_4.0.4           bookdown_0.22         ellipsis_0.3.2        spatstat.core_2.1-2   jquerylib_0.1.4       RColorBrewer_1.1-2    ggridges_0.5.3        plyr_1.8.6            rpart_4.1-15          deldir_0.2-10         pbapply_1.4-3         viridis_0.6.1         cowplot_1.1.1         zoo_1.8-9             haven_2.4.1           ggrepel_0.9.1         cluster_2.1.1         fs_1.5.0              magrittr_2.0.1        RSpectra_0.16-0       data.table_1.14.0     scattermore_0.7       openxlsx_4.2.3        DO.db_2.9             lmtest_0.9-38         reprex_2.0.0          RANN_2.6.1            fitdistrplus_1.1-5    matrixStats_0.59.0    hms_1.1.0             patchwork_1.1.1       mime_0.10             evaluate_0.14         xtable_1.8-4          readxl_1.3.1          gridExtra_2.3         compiler_4.0.4        KernSmooth_2.23-18    crayon_1.4.1          shadowtext_0.0.8      htmltools_0.5.1.1     mgcv_1.8-36           later_1.2.0           lubridate_1.7.10      DBI_1.1.1            
## [109] tweenr_1.0.2          dbplyr_2.1.1          MASS_7.3-53.1         Matrix_1.3-4          cli_2.5.0             igraph_1.2.6          pkgconfig_2.0.3       rvcheck_0.1.8         plotly_4.9.4          spatstat.sparse_2.0-0 xml2_1.3.2            bslib_0.2.5.1         rvest_1.0.0           digest_0.6.27         sctransform_0.3.2     RcppAnnoy_0.0.18      spatstat.data_2.1-0   rmarkdown_2.8         cellranger_1.1.0      leiden_0.3.8          fastmatch_1.1-0       uwot_0.1.10           shiny_1.6.0           lifecycle_1.0.0       nlme_3.1-152          jsonlite_1.7.2        limma_3.46.0          viridisLite_0.4.0     fansi_0.5.0           pillar_1.6.1          lattice_0.20-41       fastmap_1.1.0         httr_1.4.2            survival_3.2-10       GO.db_3.12.1          glue_1.4.2            remotes_2.4.0         zip_2.2.0             png_0.1-7             bit_4.0.4             ggforce_0.3.3         stringi_1.6.2         sass_0.4.0            blob_1.2.1            memoise_2.0.0         irlba_2.3.3           future.apply_1.7.0

Cytotoxic: Level 4

2021/05/13